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Abstract 

We construct an improved version of nonrelativistic QCD for use 
in lattice simulations of heavy quark physics, with the goal of reducing 
systematic errors from all sources to below 10%. We develop power 
counting rules to assess the importance of the various operators in the 
action and compute all leading order corrections required by relativ- 
ity and finite lattice spacing. We discuss radiative corrections to tree 
level coupling constants, presenting a procedure that effectively re- 
sums the largest such corrections to all orders in perturbation theory. 
Finally, we comment on the size of nonperturbative contributions to 
the coupling constants. 
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1 Introduction 



Several papers in recent years have explored the possibility of substituting a 
nonrelativistic action for the usual Dirac action in lattice QCD simulations 
of heavy quarks |l|, |], |[ . In particular, this technique has been applied to 
simulations of the ip and T families of mesons. Results from these simula- 
tions have been very encouraging. They not only capture the gross features 
of quarkonium physics but also are surprisingly accurate given the many 
approximations inherent in the simulations. Such success follows from a va- 
riety of properties unique to quarkonium states and from the algorithmic 
advantages of a nonrelativistic formulation of quark dynamics. This suggests 
the possibility of a truly reliable numerical study of quarkonium states us- 
ing computing resources that are currently available, a study in which all 
systematic errors are identified, computed, and removed to some reasonable 
level of precision. In this paper we develop the outline for such a program 
and provide some of the ingredients needed to move from existing quarko- 
nium studies to precision studies. Most of the techniques we describe are 
readily generalized for use in simulations of other bound states containing 
heavy quarks, such as D and B mesons. 

The major tools for simulating gluon dynamics in QCD all involve Monte 
Carlo methods of one sort or another. Consequently, in simulating quarko- 
nium states, we must worry about statistical as well as systematic errors. 
However, simulation data can be generated far more efficiently for quarko- 
nium states than for ordinary hadrons, and thus the dominant errors for many 
of the most important quarkonium measurements are systematic. Here we 
focus only on these. 

The nonrelativistic formulation of heavy-quark dynamics is referred to as 
nonrelativistic quantum chromodynamics or NRQCD. There are five poten- 
tially significant sources of systematic error when NRQCD is used to study 
quarkonium states: 

• Relativity. The basic action for the heavy quarks incorporates the 
physics of the Schrodinger equation. Left out are relativistic corrections 
of where v is the quark velocity. Since v"^ is approximately 0.3 

for the ip's and 0.1 for the T's, we expect errors of 10-30% due to 
relativistic corrections. Such errors can be systematically removed by 
adding new interactions to the heavy-quark action. 
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• Finite Lattice Spacing. As in all lattice simulations, there are errors 
resulting from the discretization of space and time. For our systems, 
such errors are (9(a^p^) and 0{aK), where a is the lattice spacing and p 
and K are the typical three-momentum and kinetic energy of a quark in 
the hadron. Since typically a ~ where M is the quark mass, these 
errors are roughly the same order as those due to relativity. These may 
also be systematically removed by correcting the heavy-quark action. 

• Radiative Corrections. As in the Dirac action, the basic action in 
NRQCD contains only two parameters: the QCD coupling constant g, 
and the quark mass M.0 However, each of the new interactions added 
to NRQCD to correct for relativity or finite lattice spacing introduces 
a new coupling constant. These new couplings can be treated as free 
parameters, to be adjusted by fitting to experimental data, but this 
significantly reduces the predictive power of a simulation. Alterna- 
tively one can try to compute these couplings in terms of g and M 
using perturbation theory, thereby reducing the number of parameters 
back to two. The radiative corrections to these coupling constants are 
determined by physics at distances of order the lattice spacing a, and 
therefore perturbation theory should be applicable provided a is small 
enough. However, a must not be too small. NRQCD is nonrenor- 
malizable, and as a result the radiative corrections contain power-law 
ultraviolet divergences, like g'^/aM, which make perturbation theory 
useless unless a ~ 1/M or larger. Lattice spacings currently in use, 
where (7^ ~ 1, are probably suitable for both c and b quarks, and also 
seem small enough for perturbation theory to work. Radiative cor- 
rections of order 20-40% could easily result, making these among the 
largest contributions at the order we are considering. We will see that 
they are well estimated using mean-field theory techniques. Uncalcula- 
ble corrections due to nonperturbative physics at short distances could 
be of order a few percent. 

• Finite Lattice Volume. Quarkonium hadrons are much smaller than 

^ Actually there is a third parameter, Eq, which adjusts the origin of the nonrelativistic 
energy scale to coincide with that of the relativistic energy scale. However, this parameter 
affects only the absolute masses of hadrons and has no effect upon nonrelativistic quark 
dynamics, mass splittings, wavefunctions, etc. In principle Eq may be computed, at least 
approximately, using perturbation theory. 
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ordinary hadrons. Typical lattices in use today can easily be 10 T- 
radii across, making it unlikely that the finite volume of the lattice has 
much of an effect upon simulation results. 

• Light- Quark Vacuum Polarization. Hadrons built of heavy quarks are 
affected by light quarks through vacuum polarization. Such effects are 
usually very expensive to simulate numerically due to the inefficiency 
of the algorithms when quark masses are small. However, effects from 
light-quark vacuum polarization are likely to be small for quarkonium 
states. Some evidence comes from the decay widths of excited ip and T 
mesons into D and B mesons, respectively. These decays result from 
light-quark vacuum polarization, and their widths are closely related 
to the energy shifts due to light quarks. The small size of these widths 
compared to typical mass splittings, together with model studies of this 
coupled-channel problem [Q , suggests that light-quark vacuum polar- 
ization might affect quarkonium states at only the 10% level. Further- 
more, the extrapolation from large to small light-quark masses should 
be quite smooth for low-lying quarkonium states since, unlike the p 
meson, for example, these are all effectively stable even for realistic 
light-quark masses.^ Note finally that heavy-quark vacuum polariza- 
tion has only a small effect upon quarkonium states since there is too 
little energy in such states to easily produce a heavy virtual quark and 
antiquark. For most purposes such effects can be ignored, though they 
are easily incorporated through perturbatively computed corrections to 
the gluon action. 

This list suggests that we might be able to reduce systematic errors be- 
low 10% in quarkonium simulations using NRQCD. To achieve this goal we 
must compute all of the leading-order interactions that correct the action 
for relativity and finite lattice spacing as well as the leading radiative cor- 
rections for the most important of these interactions. Other sources of error 
are probably unimportant at this level, although including some light-quark 
vacuum polarization effects (by extrapolation) is both feasible and desirable. 

A low-lying quarkonium state does decay, via the annihilation of the heavy quarks 
into gluons. The decay rate is relatively small and the decay mechanism has a negligible 
effect upon other properties of the hadron. Consequently, this decay mechanism is usually 
omitted in simulations. 
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In this paper, we present all tree-level corrections to the NRQCD action 
that are relevant to work at this level of precision. Although we have only 
just started the analysis of radiative corrections, we also present here some 
general comments concerning the nature and treatment of these corrections. 

Our analysis begins in Section ^ with the development of power-counting 
rules that allow us to determine the relative importance of various operators 
in the NRQCD action for analyses of quarkonium states. Such rules are triv- 
ial when one is considering light-heavy systems such as D oi B mesons. They 
are more complicated for heavy-heavy systems since the heavy quark is not 
a spectator in such a hadron; in particular, the quark mass strongly affects 
the hadron's size and internal dynamics. In Section ^ we compute the (9(f ^) 
corrections to the NRQCD action in the continuum for both Minkowskian 
and Euclidean space. In Section |^, we discretize the corrected NRQCD ac- 
tion, further modifying it to remove the leading errors due to the finite lattice 
spacing. In Section |^, we discuss some general properties of the radiative cor- 
rections, we give estimates for the largest of these, and we discuss the extent 
to which nonperturbative physics can affect the NRQCD couplings. Finally, 
in Section H, we summarize our results and discuss future prospects. 



NRQCD is an effective field theory which approximates ordinary relativistic 
QCD at low energies. In principle, the NRQCD action may be corrected 
to reproduce the exact results of QCD by including an infinite number of 
nonrenormalizable interactions. In practice, however, exact agreement is 
unnecessary: only a finite number of interactions are required to attain any 
desired accuracy. Some criterion is necessary to assess the importance of 
various interactions when designing the corrected action. In this section we 
categorize interactions by their effect on the energy spectrum of quarkonium 
states. 



Correction terms 6C{x) in the lagrangian density of NRQCD contribute 



2 Power-Counting Rules for NRQCD 



2.1 Building Blocks for NRQCD 




(1) 
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to the energy of a quarkonium state \n). Since operators are not normal- 
ordered in lattice simulations, there are cases where large contributions to 
5E involve contractions of two or more fields in a single term of 5C (x) . Such 
contributions have no bearing upon the relative importance of corrections 
since their effect is cancelled by shifts in the coupling constants of lower- 
dimension interactions. Thus we may ignore these contributions when esti- 
mating the importance of various terms. For the remaining contributions, 
each field and operator in 5C{x) has a characteristic magnitude determined 
by the dynamics of the hadron. 

The important dynamical scales for the nonrelativistic analysis of the 
ip and T states are the typical momenta p and kinetic energies K of their 
quarks; these are of order 1-2 GeV and 500 MeV respectively. If v is the 
quark velocity, the two scales are related to the quark mass M by 

pr^ Mv K ^ Mv'^. (2) 

These determine the characteristic magnitudes of the fundamental fields and 
operators in NRQCD and allow us to estimate the (normal-ordered) contri- 
bution of any interaction to a quarkonium energy. 

The continuum action for NRQCD is built from the fields and operators 
shown in Table We also list the estimated magnitude for each of these 
in terms of v and M. It is important to remember that it is the order in v 
rather than the dimension of an operator that determines its numerical im- 
portance. In this table, the quark field ip and antiquark field x have two spin 
components and three color components that transform in the fundamental 
representation. Combinations such as ip^tp or x'^'^ are invariant under rota- 
tions and color transformations, while ip'^crip, for example, is a color singlet 
but rotates as a vector. The gauge fields have two color indices and transform 
in the adjoint representation. For example, the chromoelectric field is 

E(x) = E"(x)T". (3) 

The generators T° are normalized such that Tr(T"T'') = 5°^/2. The combi- 
nation ip^'Eip is then a color-singlet vector field. Covariant derivatives act on 
the quark and gauge fields according to 

Dtipix) = {dt + ig4>{x))i!{x) 
Dtjj{x) = (V - igA{x))'ip{x) 
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DtE,{x) = dtE,{x)+i[g^{x),E,{x)] (4) 
BE,{x) = VE,{x)-t[gA{x),Ei{x)], 

where and A are the scalar and vector gauge potentials, and g is the 
coupling constant. The last definitions insure that the chain rule works 
properly for covariant derivatives. For example, 

D{EiIj) = (DE)^ + E{D^Ij). (5) 

The nonabelian electric and magnetic fields can be defined in terms of these 
covariant derivatives: 

[D,,Dj]ij = -ige^^^kB^^i^. (6) 

The estimated sizes given in Table follow from the properties of quarko- 
nium states. For example, the number operator for (heavy) quarks, 

d^xi)\x)il){x), (7) 

has an expectation value very near 1 for a quarkonium meson. Since the 
quark in the meson is localized within a region Ax ~ we estimate 

pi 

and therefore 

(9) 

that is, %Ij Similarly, the operator for kinetic energy, 

d^xi,\x)^i,{x), (10) 

has an expectation value of K by definition, and so the spatial covariant 
derivative acting on a quark field is of order 

T> ^ {2M Kf'^ ^ (11) 

as expected. 
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Field equations can also be used to relate estimates for different operators. 
For example, the lowest-order field equation for the quark field, 

implies that 

K, (13) 

again as expected. If we specialize to Coulomb gauge, the gauge most natural 
for nonrelativistic systems, this equation becomes 



2M 



+ — 1 V~0, (14) 



where we are neglecting the vector potential, which is small in this gauge, as 
discussed below. The potential energy that balances the kinetic energy and 
produces a bound system enters through the operator g(f) and we expect that 

g(j>{x) ~ K (15) 

in Coulomb gauge. 

This result can be checked against the field equation for (p (again neglect- 
ing the vector potential), 

V''g(t)(x)^-g^i:\x)ij{x), (16) 

which implies that 

g(l){x) o9^P^^9^P- (17) 

This is consistent with the previous estimate if the effective low-energy cou- 
pling strength is 

Ois ~ 9'^ ~ f ■ (18) 

That the Coulomb-gauge vector potential in a quarkonium state is typi- 
cally smaller than the scalar potential may be inferred from the field equation 
for the vector potential, 

{d^ - V^) = ^i^^Vi: + g^Vg^ + • • • . (19) 
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Using the estimates given above, we see that 

gA{x) r^^^{^p^+ pK'^ ^ vK (20) 

which is smaller than the scalar potential by a factor of the quark velocity. 
Estimates for the nonabelian electric and magnetic fields follow immediately 
from 

g'E — —'Vg(j) + ■ ■ ■ ~ pK 

gB = V X gA + -- - ^ . (21) 

In quarkonium, as expected for a nonrelativistic system, magnetic fields are 
smaller than electric fields by a factor of v. 

These estimates are valid perturbatively, as in the analysis of positron- 
ium with NRQED, but are also consistent nonperturbatively. The essential 
ingredient in their derivations is the nonrelativistic nature of quark dynamics. 



2.2 NRQCD interactions in the continuum 

Equipped with the power-counting rules of the last section, we can enumerate 
terms for the quark action in NRQCD that are potentially important for 
quarkonium physics. The leading terms are those of the Schrodinger theory: 




Correction terms due to relativity and finite lattice spacing must respect the 
symmetries of the theory: gauge invariance, parity, rotational symmetry, uni- 
tarity, and so on. For example, the interaction ip^E-cril), corresponding to an 
intrinsic electric dipole moment, is not allowed because it is odd under parity, 
while ip'^'B ■ (Tip is even under parity and so is allowed. Charge-conjugation 
invariance requires that the total action be invariant under the interchange 
ip X- Thus the antiquark action can be obtained directly from the quark 
action. 

Correction terms must also be local. This, combined with the fact that 
our theory needs only to be accurate to 10%, severely restricts the number 
of interactions. For example, the interaction ip^'B'^ip/M^ is consistent with 
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required symmetries and so must appear in the NRQCD action. But it is 
unnecessary for our calculations because it is suppressed by (~ 10^^ for 
T's) relative to the leading terms in the action. 

Ignoring spin splittings for the moment, only corrections that are sup- 
pressed by f ^ relative to the leading terms are needed to achieve accuracy to 
10%, at least for T's. The only such terms bilinear in the quark field are the 
order v"^ corrections 

(^>Cbilinear = ^1 Ip^B^lp 

+ C2^^t(D.E-E-D)^ 

+ C3^V'^o--(DxE-ExD)V' (23) 

In this and subsequent equations, D acts on all fields to its right. The dimen- 
sionless coefficients ci . . . C4 are functions of the running coupling constant 
as{TT/a) and aM. 

Note that terms involving time derivatives of the quark field, such as 
ip^D^tlj, are not included. Such terms greatly complicate the numerical eval- 
uation of quark propagators. We may avoid them here by suitably redefining 
the quark field so that factors of iDt are in effect replaced by factors of 
— D^/2M, in accordance with the field equation for ip, 

iDM^)^^^{x). (24) 

These transformations do not affect masses, on-shell scattering amplitudes, 
or other physical quantities. 

In addition to the bilinear terms, there are four-fermion contact interac- 
tions involving a quark and antiquark: 

+ ds^V'Vx-xW- (25) 

These appear to be down by only a single power of v. However, similar in- 
teractions do not occur in relativistic continuum QCD, and therefore such 
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terms can occur here only to one-loop order and beyond. Thus the coeffi- 
cients di and d2 are both of order aKir/a), making these contact interactions 
significantly less important than the bilinear interactions considered above. 
Contact interactions also occur between pairs of quarks, but these can only 
affect heavy-quark baryons. Four-fermion operators can also couple to col- 
ored states: 

scholar ^ ds^Y^^^T^XX^Ty 

+ d.^Y.^^T'^^x-X^T'^fTi;, (26) 

where d^ and d^ are of order a^(7r/a). Color-singlet mesons are affected by 
these interactions, since the meson can become colored by emitting a virtual 
gluon. Gluon emission however is suppressed by f ^, so these interactions are 
again not as important as the bilinear terms in Eq. (pSj). 

Nontrivial spin dependence first appears in the bilinear correction terms 
of Eq. (^), so that spin splittings in quarkonium meson families are sup- 
pressed by f^. A 10% determination of these splittings therefore requires the 
retention of spin- dependent interactions suppressed by v^. The additional 
quark-bilinear terms required to this order are of three sorts: permutations 
of the operators ip^crip, B, and two D's; permutations of the operators ip'^crip, 
D X E, and two D's; and the operator ip'^cr ■ E x E?/', which is nonzero in 
nonabelian theories. Of these, only three interactions occur in our treatment: 

+ /2^V^t{D2,o--(DxE-ExD)}V' (27) 
+ /a^V^V-ExE^. 

The power-counting rules indicate that spin splittings in quarkonium sys- 
tems should be smaller than splittings between radial and orbital excitations. 
This is a familiar feature of the experimental data. The ip'-ip and T'-T mass 
splittings are both around 600 MeV. Our analysis indicates that spin split- 
tings should be smaller than this by a factor of f ^; that is, approximately 
180 MeV for t/^'s and 60 MeV for T's. In fact experiments show that the 
p-state hyperfine splitting, x(2"'"+)-x(0^"''), is 141 MeV for the if) family, and 
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34 MeV for T's. The sizes of these sphttings are consistent with our expec- 
tations, and show clear evidence of appropriate scahng with v"^. The s-state 
hyperfine sphtting between the if) and ?7c of about 117 MeV is also consistent. 
These data reinforce our confidence in the power-counting rules, and in the 
utility of a nonrelativistic framework for studying these mesons. 

One final source of corrections comes from vacuum polarization of heavy 
quarks. These modify the gluon action at order as{M) f ^ and so can probably 
be ignored at the 10% level. Such corrections are essentially relativistic 
in character and should not be analyzed using NRQCD. They are easily 
computed using continuum perturbation theory. 

3 Relativistic Corrections 

The corrections to the continuum NRQCD action required by relativity were 
enumerated in Section |[ Our task now is to find values for the coupling 
constants multiplying these operators such that the predictions of NRQCD 
agree with those of ordinary QCD through order f^. Since the coupling 
constants have perturbative expansions, we will determine them by matching 
perturbative results in the two theories. Here we work only to lowest order, or 
tree level, in perturbation theory. Higher order corrections will be discussed 
in Section ^. 

3.1 Kinetic Terms 

The simplest correction follows immediately from the formula for the rela- 
tivistic energy of a noninteracting quark. 




(28) 



This implies a correction term, appropriately gauged, of the form 



5Ck 



1 



(29) 



8M3 



which fixes Ci = 1/8 in Eq. (pS]). 
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3.2 Electric Couplings 



Most of the corrections from Section ^ involving the chromoelectric field are 
linear in gE. We compute these terms by examining the order g amplitude Te 
for scattering a quark off a static electric field in QCD, 



Te{p, q) = u{q)-f° ^</)(q - p) u{p) 



(30) 



with the scalar potential. We are interested in matching this result at small 
V to that obtained in NRQCD, and so must expand it in powers of q/M and 
p/M. The Dirac spinor, normalized nonrelativistically with u'^u = 1, is 



u{p) 



2Er, / 



• P 

Ep + M 



(31) 



with Ep = -y/p^ + M'^ and ip a two-component spinor. Substituting this 
expression in Eq. (|30D we obtain 



TE{p,q) 



iEp + M){Eg + M) 



1 + 



AEpE, 

p • q + icr ■ q X p 

{E, + M){Ep + M) 



g(j){q-p)'^ 



(32) 



= ^i?(p,q) + VB(p,q) , 

where Se denotes the scalar part of the amplitude, independent of cr, and 
Ve is the spin-dependent term. Expanding both Se and Ve to second order 
in p/M and q/M, 



'^i?(p,q) 

VE(p,q) 



1 - 



IP - qj^ 

3i 



ip^ g(f){q_ - p) ^ 



(33) 



4M2 32M4 



P + q)U'^cr-qxp g(j){c[ - p) i> 



By computing the same process in NRQCD and comparing with the QCD 
result at low momentum, it is possible to determine which operators must 
appear in the NRQCD lagrangian and to fix their coefficients. The first term 
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in the lowest order NRQCD action, Eq. (p2D, predicts a scattering amplitude 
ip"^ g(j){q — p)iIj . This corresponds to the first term in Se- The remaining 
terms in Se and Ve can only be included by adding the new interactions 

+ ^^^(^T-DxE-o--ExD)^ (34) 
+ ^^V^t{D2,^T-DxE-o--ExD}^. 

to the NRQCD action. The coefficients C2 and C3 in Eq. ( p3D are therefore 
1/8, while /a in Eq. (0) is 3/64. 

3.3 Magnetic Couplings 

Similarly, we determine the operators linear in B by calculating the amplitude 
for the scattering of a quark off a static vector potential A(x). The time 
independence of A(x) guarantees that there is no electric component. It also 
enforces the conservation of the quark energy. 
The amplitude is 

^i?(p,q) = -M(q)7-^A(q-p)M(p), (35) 

where = q^. Expanding in p/M and employing energy conservation. 



2M \ 2A/P J 
X ^/'"^ [(p + q) • A + icr ■ A X (p - q)] (36) 



= 5'ij(p, q) + VB(p,q) • 

The spin- independent term Sb{p, q) arises in NRQCD from the terms linear 
in gA which appear in the kinetic energy part of the action, Eqs. (|22| ) and 
( ^91) . The spin-dependent term Vb must be generated by the new interactions 

+ ^^HD^^■B}V., (37) 
so that C4 = 1/2 and /i = 1/8 in Eq. (|3|) and Eq. (gT]) . 
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3.4 Bilinears in F^jy 

The interaction ip'^ cr ■ E x E ip , which contributes to spin sphttings, is the 
only bihnear in F^^ that contributes to the accuracy desired. To fix the 
coefficient of this operator (/a in Eq. (^)), we calculate the amplitude for 
double scattering of a quark off an external static electric ffeld. In QCD, it 
is 

TEEiPu P2) = / (27r)=^53(qi + + Pi - P2) (38) 

X M(p2)7°#(q2) 7;^— — T^-— #(qi) 7°w(pi) • 

{pi + ii) - M + ie 

This same process must be computed in NRQCD, and compared to an ex- 
pansion of Eq. ( P^D at small p/M. This comparison is simpliffed by using 
the identity 



k^-E k^ + E 



(39) 



for the fermion propagator, to break Tee into two terms before expanding. 
NRQCD, with the terms computed thus far, reproduces the first term con- 
taining u u/ (k^ — E) to the required order in v, using vertices derived from the 
interactions of Eq. (|3^). These were designed to match the matrix element 
of Eq. (PPD, M7°(y(0M, and the first term of Eq. (^) contributes the same 
matrix elements to each of the vertices of the QCD diagram. The energy 
denominator 

1 1 _ 1 / kM 1 

k^-E^ k^- kV2M ~ ^0 _ k2/2M [sM^J fco - kV2M ' ^ ' 

with the kinetic energy k^ = k^ — M, is reproduced in NRQCD by its non- 
relativistic propagator and the kinetic term correction of Section 

The second part of Tee, when expanded, introduces a term of an entirely 
different character. To lowest order in v, for low-momentum external gluons, 
the antiquark propagator is just (2M)~^, so this part behaves as a local 
seagull interaction. This contributes to Eq. ( p8|) 

X ^^^^ ^0(q2) [qi ■ q2 - io- ■ qi X q2]^0(qi) ip ■ 
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In order to obtain the same result with NRQCD, we must include the new 
interaction ^ 

^Cee = -^^^ [E ■ E + ^o- ■ E X EjV' . (42) 

To the order to which we are working only the second term, which contributes 
to spin splittings at order f^, must be retained. So /s = —1/8 Eq. (|27|), 
and we have determined all relevant tree-level coefficients of the NRQCD 
lagrangian. 



3.5 The NRQCD Continuum Lagrangian 

The lagrangian for continuum NRQCD in Minkowski space is 

Cnrqcd = tipWiifj + ifj^ —t/j 

+ SCk + SCe + SCb + SCee, (43) 



where the relativistic corrections are defined in Eqs. (pOf), (|34D , (p^) and 
respectively. We have not included in the lagrangian a mass term Mip'^ip be- 
cause in a nonrelativistic framework it only fixes the zero-point energy and 
has no effect on mass splittings, wavef unctions, and so on. As an independent 
check, we performed the systematic but lengthy calculation of the lagrangian 
in Eq. (^) using the Foldy-Wouthuysen-Tani transformation |^ . This trans- 
formation is not directly applicable beyond tree level. However, perturbative 
corrections to the various coefficients may be computed by comparing am- 



plitudes, as in Sections 2]2| - 13^ , but including loop corrections. 

Finally, we note that for lattice simulations we will need the Euclidean 
action. This is obtained from the Minkowski theory defined by Eq. (^) 
by keeping track of how the three-dimensional vectors and scalars are de- 
fined in terms of Lorentz-covariant quantities. Under Wick rotation the zero 
component of contravariant vectors rotates as the time coordinate, 

^(M) = ~'^^\e) 1 (44) 
while the zero component of covariant vectors rotates as the time derivative, 

= ) . (45) 
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Spatial components of all four- vectors are unchanged. Wick rotation changes 
the sign of the (0, 0) component of the metric tensor, so that the usual 
Minkowski metric rj^^i, = diag(l, — 1, — 1, — 1) rotates to the negative of the 
identity matrix. One then has to extract this minus sign from scalar prod- 
ucts in order to work with the usual positive definite euclidean metric. Our 
definitions for the basic operators appearing in Eq. (^31) are 

4>{M) = -i4>{E) (46) 
while other fields are unaffected. 



4 Lattice NRQCD 
4.1 Leading Order 

NRQCD is readily reformulated on a discrete space-time lattice. As usual, 
quark fields ip{x) are defined at the nodes of the lattice, while gauge fields are 
represented by unitary matrices Ux,^ defined on the links joining neighboring 
nodes. Covariant derivatives are replaced by forward, backward or centered 
differences, 

A{±) = i(AW+A(-)), (47) 
depending upon the details of the interaction, while the laplacian becomes 
A(2)^5^aWa;-)=EaS-)aW. (48) 

i i 

The field F^y{x) is efficiently represented by cloverleaf operators defined at 
the nodes of the lattice , 

9Ft^K^) = -^, E 2:[f/p(.,,.)] , (49) 
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where the sum is over all plaquettes P in the (/i, v) plane containing the site x. 
Up is the product of link matrices on P, counterclockwise about fi x D, as 
depicted in Fig. 1. The operator X is defined as 

X[M] = ^ _ 1 ini(Tr M) . (50) 

2z 3 

Covariant derivatives of F^u are represented by the differences 



x] 



aA(-)FW(a;) ^ ^^(x) - t/t,^,,FW(x - ap)f/._.^,, . (51) 

All terms in the NRQCD action are built with these elements. 

Focusing on the leading-order terms, we may define a series of lattice 
actions for NRQCD, 

= a'Y.^\x),p{x) 

X 

- a'j:^\x + ai)(l-^y Ul(l-^)\ix), (52) 

where n is a positive integer and Hq is the kinetic energy operator 

Ho = — . (53) 

2M ^ ^ 

The quark Green function for such a theory satisfies a simple evolution equa- 
tion, 

^(x, t + a) = [l-^j Ul [1 - ^ j G(x, t) + 5^,oSt+a,o , (54) 

where G(x, t) vanishes for t <0. This evolution equation differs from others 
that have been suggested. It is symmetric with respect to time reversal, un- 
like that of Also, as we will see, it leads to a wavefunction renormalization 
much smaller than that of 0. 

The parameter n was introduced to prevent instabilities at large momenta 
due to the kinetic energy operator [Q. The Green function defined by Eq. (p^) 
blows up if max{aHo/2n) > 2, and so n must be large enough to avoid this. 
Neglecting gluons, this implies that n > 3/(2Ma) is required for stability; 
gluonic effects reduce the 3/2 to something closer to 1.15 for /3 ~ 6 p. Thus, 
at /? = 6, n = 1 suffices for b quarks, while n = 2 is needed for c quarks. 
Larger n's may be required as /5 is increased and a decreases. 
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4.2 Lattice Spacing Errors 

The finite spacing of tlie lattice introduces new systematic errors into NRQCD. 
Errors associated witli tlie temporal spacing at are typically of order at{K), 
while those due to the spatial lattice are of order a^(p^). The two tend to 
be roughly comparable since usually ax = at ^ 1/M. Both must be removed 
if we are to achieve high precision. 

Errors due to the finite lattice spacing are removed by adding new in- 
teractions to the lagrangian, just as 0(f ^) errors were removed. As in that 
case, we could compute the necessary corrections by matching perturbative 
scattering amplitudes in lattice NRQCD with those in continuum QCD. How- 
ever, at tree level, it is simpler to correct the individual components from 
which the lattice action is built (A*^^\ E, B. . . ) so that they more accurately 
reproduce the effects of their continuum analogues. This is the approach we 
will follow here. 

To compare the effects of lattice operators with those of continuum op- 
erators we need some way to relate fields on the lattice to those in the con- 
tinuum. For the quarks, we simply identify the lattice field at a site with 
the continuum field at the same space-time point. For the gluons, gauge 
invariance requires that the link operator be related to the continuum gauge 
field through the path-ordered exponential of a line integral. 



where the scalar product is taken with the positive euclidean metric. The 
simplest choice of path for the integral is a straight line joining x to a; + a/i; 
other choices lead to the same final results but give more complicated 0{a) 
corrections. With this mapping from lattice variables to continuum variables, 
we can now correct each of the lattice operators. 

4.2.1 Spatial Derivatives 

Given relation (|55|) , the action of the lattice difference operators on contin- 
uum fields is specified by 




A ■ dy 



(55) 




(56) 
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These are just gauge- covariant extensions of the obvious g = relations. By 
combining these expansions we obtain an improved difference operator that 
reproduces the effects of Di through order a"^: 



(57) 



Similarly we may define an improved lattice laplacian that is also accurate 
to order a^: 



12 ^ L * 



,(-) 



(58) 



4.2.2 Temporal Derivative 

Temporal and spatial derivatives enter NRQCD differently because the the- 
ory is nonrelativistic. One consequence is that quark propagation is governed 
by a Schrodinger equation with a single time derivative. The calculation of 
quark Green functions is therefore an initial value problem, which is much 
less costly to solve numerically than the boundary value problem dictated by 
the Dirac equation. If we improve the time derivative operator as we did for 
spatial derivatives, this simplification is lost, as we would need to introduce 
higher order time derivatives. 

To find an alternative, we examine the evolution equation for the quark 
Green function, Eq. (|54D . Neglecting the gauge field for the moment, this 
equation may be written (for t > 0) 

G(x,t + a) = (l-^) G(x,t) 

= e-'^^^ffG(x, t) , (59) 

where the effective hamiltonian 



a \ 2n 

The order a term in the effective hamiltonian cancels if we replace Hq by 

Ho^Ho- ^Hl . (61) 
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Making this replacement in the full evolution equation (Eq. (|5^)) removes the 
leading error due to our lattice approximation of the temporal derivative.^ 
The resulting theory is equivalent to what would have been obtained by 
improving the temporal derivative along the lines presented in the previous 
section; the two theories are related by a redefinition of the quark field, and 
so give the same results for energy levels, on-shell scattering amplitudes and 
other physical quantities. 



4.2.3 Electric and Magnetic Fields 

The lattice cloverleaf field F'^J{x) (Eq. (|49|)) is equal to the continuum 
field F^y{x) up to corrections of order a?. Our power-counting analysis (Sec- 
tion ^ implies that only corrections linear in F^^, are important. Such cor- 
rections are the same in an abelian theory as they are in QCD, so we simplify 
the analysis by focusing on the abelian case. 

In the abelian theory, where the operator X just takes the imaginary part, 
the cloverleaf field is 



F^^Jix) = -Im(l-^i A-dy + 



' \ 4 /(2X2) 



/ {d,A, - d,A^) dy^ A dy'' + 0{a') (62) 

J(2x2) 



9_ 

4 7(2x2) 



= a'gF,y{x) + -{dl + dl)gF,y{x) + 0{a''), 

where the surface integral is over the (2 x 2) plaquette in the (/x, z/) plane 
centered at x. The term comes from the second-order term in the Taylor 
expansion of F^y around point x. The nonabelian generalization of this result 
is obtained by replacing derivatives with covariant derivatives. 
By subtracting the lattice version of the 0{a^) error, 

gFj^Kx) = gFlfJix) - J [aWa(-) + A(+)A(-)] gF!fJ{x) , (63) 

we have a definition accurate to O(a^). Expanding the lattice derivatives as 
in Eq. (^Tj), we obtain finally 

gF^^^Jix) = ^gF^'^Jix)-^[u,,,gF;fJix + ai,)Ul^ 



^Only the kinetic part of the hamihonian needs fixing here. The gauge-potential part 
is automaticahy exponentiated since the gauge fields enter through the link variables Ux.^- 
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+ Ul_af,,^9Fl^'i}{x-aix)a,^ 



(64) 



We have verified that F^l! reproduces the continuum F^^, to C(a^) directly 
in the nonabelian theory using techniques described in the Appendix. 



4.2.4 The Gauge Field Action 

The gluon action has lattice spacing errors as large as those of the nonrela- 
tivistic quark action. A great deal of sophisticated work has been done on 
improving the gauge field action, dealing mainly with one-loop contributions 
[0]. For us, however, a straightforward improvement of the classical action is 
all that is required. In particular, corrections involving nonplanar loops are 
not needed. This greatly simplifies our analysis. 

Power counting indicates that the nonabelian and abelian corrections have 
the same form. By Stokes' theorem, a single abelian plaquette centered at 
point X is related to the continuum field tensor by 



Up{x) = exp 



-zga%^{x) - ig-{dl + dl)F,^{x) + 0{a') 



(65) 



and the usual lattice lagrangian is 

^ Re(f/p(a:)-l) = -lF^2^(a;)-^F,.(a:)(9; + 9^)F,.(x) + (9(a^). (66) 



The nonabelian generalization is straightforward; the key result is that the 
corrections are of order a? . 

A simple scaling argument allows us to correct the nonabelian action. 
The usual lagrangian involves a sum over all four (1x1) plaquettes at a 
site X, 

1 



C-g{x) 



Tr[7^(f/p(lxl) 

P{lxl) 



where 



7^(M) = 



M + Mt 



1 . 



From the discussion above. 



(67) 



(68) 



(69) 
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where is the continuum lagrangian. An acceptable lagrangian may also 
be constructed from (2 x 2) plaquettes, and its a-dependence follows from 
Eqs. (^rp and with the replacement a — > 2a: 

(f^P{2x2)) 
""i* P(2x2) 

= ^g'"* + + C(a^) • (70) 

By combining these two lagrangians, we obtain a gluon lagrangian that is 
accurate to order a'^, 

+ O(a^) . (71) 

Again, this result may be confirmed directly in the nonabelian theory by 
means of the techniques described in the Appendix. 
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4.3 A Corrected Evolution Equation 



The leading-order evolution equation, Eq. (|5J), for the quark Green function 
is readily generalized to include relativistic and finite a corrections. We 
choose the full evolution equation to be 

G(x,t + a) = 

(-#)"(-1^)<'(-1^) 

for t > 0. Here we have corrected the leading kinetic energy operator for 
finite lattice spacing errors as in Eqs. (|58|) and ( |6TD so that 

~ _ A(2) a (A(2))2 

^° = "2M "i^^^M^- ^^^^ 

Further, a straightforward transcription of the relativistic corrections into 
lattice variables gives, for the spin-independent corrections, 

^ ^(^^^^-E-E-AW), (75) 
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and, for the spin-dependent terms, 



^^spin = -^^-(AWxE-Ex AW) 

^ cr-B (76) 



2M 

^i^spin,.^ = -^{A(2),^.B} 



{A(2),o--(AWxE-Ex aW)} 



64M4 

- ^cr-ExE. (77) 

Here E and B are defined in terms of the cloverleaf field, Eqs. (^) and (p^), 

E^(a:) = fI^\x) (78) 
i?Xx) = '^U,kFjl\x). (79) 

The difference operators in these terms act on all the fields, E, B, and G(x, t), 
to their right. However, the D ■ E — E ■ D operator may be simplified since 
the net effect is to differentiate only the E field. Using Eq. (|5l|) , it can be 
written as 

^^^^ ■ E = ^ E [ U.,E\x + ai)Ul - Ui_,,^,E\x - ai)U,.ai, ] • (80) 

Finite lattice spacing corrections are included in 5-ffspin by using E and B 
of Eq. (H) and A(±) of Eq. in place of E, B and A(±). These are 
unnecessary in the other operators, which are already sufficiently accurate. 

Certain correction terms may be unimportant in some calculations. The 
Oiy"^) corrections for spin-independent quantities are all included in 6Hk,v^, 
SHy2, and ^i^spin; the term 5ifspin,t,2 enters only at order f^. Spin splittings, 
however, are themselves O(f^), and so a comparably accurate measurement 
of such a splitting would require all of the corrections. For simulations of 
light-heavy mesons, such as 5's or D's, probably only the magnetic part of 
SHgpin is important. All other terms are suppressed by at least an additional 
power of the heavy quark velocity. 

For reasons of efficiency we have separated the relativistic corrections SH 
from the leading order kinetic energy Hq in Eq. (|72D. This is possible because 
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terms proportional to E or B do not cause instabilities at high momenta; it 
is only the kinetic terms that cause problems. The relativistic correction to 
the kinetic energy 511^^2 is a special case. It probably should not be in- 
cluded in the evolution equation at all, but rather its contribution computed 
in first-order perturbation theory. The problem with this term is that it is 
negative and large for large momenta, eventually overwhelming the leading 
kinetic term Hq and destabilizing the theory. This problem is not particular 
to the lattice version of NRQCD: it is a complication due to the nonrenor- 
malizability of the theory. In fact, the continuum hamiltonian 

(p+^ (p + gAf 
^ - 8M^ ^^^^ 

is not bounded below and will lead to problems if p is allowed to become 
too large. This problem does not arise in first-order perturbation theory, 
so that a perturbative analysis of this correction is possible. Note also that 
in Coulomb gauge the correction can be approximated by — p*^/8M^. Since 
this is independent of the gauge field, its perturbative contribution to mass 
splittings can be reliably determined given only the quark wavef unctions, and 
these are easily computed in a simulation. Alternatively, one could spread 
A'-^-* for this term over, say, five lattice spacings in each direction instead of 
three. This would reduce the largest effective momentum by a factor of two, 
thereby reducing the high-momentum value of SHk,v^ by a factor of 16 while 
leaving the low-momentum behavior largely unchanged. 

Many of the same problems could arise from the finite lattice spacing 
corrections in Hq. However it seems that these do not cause problems for 
relevant values of Ma. Indeed these terms seem likely to reduce the maximum 
value of Hq. 

A more radical way of dealing with instabilities at high momenta might 
be to choose Coulomb gauge, where the vector potential is small, and then 
replace 

(l - - exp (V^^T](^- m)] (82) 

in the evolution equation. Here p^ is the quark momentum, which is easily 
introduced into the simulation through the use of Fast Fourier Transforms. 
In this way the exact relativistic energy is incorporated into the simulation 
in a way that avoids instabilities. Unfortunately this approach greatly com- 
plicates the perturbative analysis of the theory, not least because it breaks 
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gauge invariance. Nevertheless such problems are probably not insurmount- 
able, and the approach may someday be useful. 

Finally, note that four-fermion interactions as in Eq. (^^) may be incorpo- 
rated naturally into the evolution equation by replacing them with new ones 
which are quadratic in fermion fields. This is accomplished by introducing 
an auxiliary field rj, a complex matrix field which is 2 x 2 in spin and 3x3 
in color indices. For example, if t] appears in the lagrangian as 

= aijj^rjijj + bx^v'^X + cTirj^r] (83) 

without a kinetic term, this is equivalent to the interaction 

-^contact = -^^'^XXV (84) 

for the fermions. This may be seen by performing explicitly the gaussian 
path integral over t], or from the equations of motion 

V = --XX^ r{^ = --i)^ljK (85) 
c c 

In a simulation this amounts to generating and summing over a random 
gaussian distribution for rj at each x and t, normalized such that 

< rjijVki >= < riij7]l^ >= SikSji , (86) 

and propagating the quark in this field according to the interactions in C^i- 

5 Radiative Corrections 

In previous sections we derived tree-level values for the couplings of NRQCD. 
The couplings are modified by radiative corrections that are dominated by 
momenta of order ir/a or larger. Since n/a is typically several GeV in sim- 
ulations, we should be able to compute these radiative corrections using 
weak-coupling perturbation theory. However, lattice perturbation theory is 
notorious for its poor convergence, which is much worse than continuum 
QCD at comparable momenta. The large corrections, which come from tad- 
pole and related diagrams and spoil lattice perturbation theory, apparently 
result from the nonlinear connection between lattice link variables and the 
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continuum gauge potential. In the worst cases such tadpole corrections are 
almost as large as tree-level contributions, and some sort of nonperturba- 
tive treatment becomes necessary. Fortunately, tadpole contributions are 
well described by a mean-field approximation [§]. Such an approach allows 
us to compute most of a radiative correction nonperturbatively in terms of 
quantities that are easily measured in simulations. In this section we first 
outline a mean-field analysis for lowest-order NRQCD, giving estimates for 
all renormalization effects in that theory. We also compare the mean-field 
theory results with exact results to first order in perturbation theory. Next 
we discuss a simple modification of lattice NRQCD that effectively removes 
all tadpole contributions. Finally we comment briefiy on the significance of 
nonperturbative contributions to the couplings. 

5.1 Mean Field Theory and Lowest-Order NRQCD 

A simple way of tracing the effects of tadpoles is to replace the link operators 
Ux,^ in the NRQCD propagator by a number uq representing the vacuum 
expectation value of Ux^p.- One gauge-invariant definition of uq is in terms of 
the plaquette operator 



Other definitions are possible, but these give similar results and will not be 
considered here. With this definition, uq is easily measured in a simulation: 
for example, uq = 1 — 0.122 at /5 = 6. The measured uq includes tadpoles 
to all orders of perturbation theory, and perhaps also some nonperturbative 
effects. 

Having a mean value for Ux^^, we may proceed with a mean-field anal- 
ysis of the NRQCD propagator. In this approximation, the kinetic-energy 
operator is 




(87) 



which becomes perturbatively 



Mo = 1 - 0.083^2 . 



(88) 




(89) 



so that 



Ho^ ho + Uq 




(90) 
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for a quark with a low momentum p, where 



/lo = 3 (1 - no) /Ma^ 



(91) 



Substituing this into Eq. (|5^) for the Green function and replacing the link 
variables by uq's, we obtain finally a mean field approximation to the quark 
propagator 



t/a 



1 - 



~2n 



auo p 

2n 2M 



2nt/a 



(92) 



In the low- momentum limit, we expect the continuum propagator to have 
the form 



Matching these, we obtain 



-t I Ef + 



E, 



mf 



In Mo (1 — a/io/2n 



2n 



for the zero-point energy induced by tadpoles. 



7mf 



Un 



[1 — ahQ/2n) 



for the mass renormalization, and 



lymi 



(93) 



(94) 



(95) 



(96) 



for the quark field renormalization. At /3 = 6 with Ma = 2 (roughly the b 
quark mass) and n = 1, these renormalization parameters are 



Ef 



Z 



mf 
M 



0.9 GeV 

1.04 

1. 



(97) 



None of the renormalizations is particularly large for leading-order NRQCD, 
at least for (5 near 6. This is despite the theory's nonrenormalizability and 
consequent power-law divergent renormalizations: for example, a/iQ oc 1/a. 
These divergences cause problems as the lattice spacing a is reduced, but 
there will be little need to reduce a once our finite-a corrections have been 
incorporated. 
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Renormalization parameters like Eq, Zm, and are needed to interpret 
simulation results. For example, the total mass of a meson with NRQCD 
energy is 

M„ = 2 {ZmM - Eo) + En . (98) 

Knowing Zm and Eq from mean-field theory, perturbation theory, or both, 
and En from simulations, we can use this expression to tune the NRQCD 
bare quark mass. The wave function renormalization is important when 
designing quark operators to describe such things as radiative transitions or 
decays through quark annihilation. For example, the continuum operator 
X^'ip, describing quark-antiquark annihilation, is well modeled by the lattice 
operator x^i^/Z^. 

We can calibrate the reliability of mean-field theory by comparing mean- 
field results with exact calculations to first order in perturbation theory. 
Such calculations have been performed by Davies and Thacker , who have 
analyzed the NRQCD propagator through one-loop order. They defined the 
propagator by means of the evolution equation 

G(x,t) = (t<0), (99) 

which differs slightly from our equation. The energy shift and mass renor- 
malization for this equation are identical to ours if n is replaced by n/2 
in our equations. The wavefunction renormalization however is significantly 
different: in mean- field theory, one finds Z^^ = (1 — a/io/"^)"" for the Davies- 
Thacker equation. In Table 2 we list the 0{g'^) coefficients for £"0, Zm-, and 
Z^. We give both the exact result, as computed by Davies and Thacker us- 
ing perturbation theory, and the mean-field estimate, obtained by replacing 
uq with its perturbative expansion, Eq. (^81) , in mean- field expressions for 
the renormalization parameters. Since we expect nontadpole contributions 
of one or two times ag/vr, or O.OSg'^ to 0.05(7^, the agreement is excellent. A 
nontadpole contribution of this size is only 5-10% of the tree-level contri- 
bution at /3 = 6. Mean-field theory does seem to account well for radiative 
corrections when they are large. 

When perturbative results are known, they can be used to correct the 
mean-field prediction. For example, the zero-point energy Eq for our propa- 
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gator at aM = 1.5 can be written as 



-In 



Mo (1 - aho/2y 



+ 0.07 



(100) 



The correction 0.07^^^ is the difference between the perturbative and mean- 
field result in Table 2. The factor 1/mq in the correction accounts for tadpole 
induced renormalization of the bare coupling g'^ = 6/j3, also suggested by 
mean-field theory ||^. With this expression, all that is needed to complete 
the calculation is a measurement of the expectation value of the trace of the 
plaquette to determine uo- 



5.2 Removing Tadpole Contributions 

Our mean-field analysis of lowest-order NRQCD is readily extended to in- 
clude the various corrections due to relativity and finite lattice spacing. There 
is however a simpler way of dealing with the tadpoles. By dividing every link 
matrix by uq before computing quark propagators, all tadpole contributions 
are automatically removed from the simulation. Using Ux,^/uq in place of 
Ux,^i everywhere, the tree-level couplings we have computed should be accu- 
rate to within corrections of order as/i^'-, that is, to 5-10% at /3 = 6.0. 

The mean-field corrections introduced by our simple procedure are not 
always small, particularly for operators that involve the cloverleaf definitions 
of electric or magnetic fields. These operators are modified by a factor Mq ^: 
for example, B B/mq. At /9 = 6.0, Uq'^ is 1.7, almost doubling the E 
and B fields. Leaving out such factors leads to dramatic underestimates of 
quantities, such as spin splittings, that involve one or more powers of E or 
B. 



5.3 Nonperturbative Effects 

As is true of all perturbative calculations, our analysis of NRQCD couplings 
omits nonperturbative contributions. A priori, we know little about the na- 
ture of such contributions. The one thing we may assert with some confidence 
is that nonperturbative contributions are smaller than perturbative contri- 
butions for momenta vr/a larger than a couple of GeV; that is, for (3 > 5.7. 
This is the fundamental assumption underlying all applications of perturba- 
tive QCD, whether in the continuum or on the lattice. As we have argued, 
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perturbative corrections to our couplings are probably less than 10-15% for 
reasonable /3's, and so it seems likely that nonperturbative effects are no 
larger than a few percent. Since most physical results are not extremely sen- 
sitive to the values of the couplings, the nonperturbative effects are probably 
safely neglected. 

Our mean-field analysis supports this view by providing a toy model for 
nonperturbative effects. The plaquette operator plays a key role in this anal- 
ysis since it determines the mean field mq. The plaquette is one of the few 
operators whose nonperturbative behavior is somewhat understood. The ex- 
pectation value of a large Wilson loop with area A contains a factor exp{—aA) 
due to nonperturbative confinement. Here a is about 0.18 GeV^, based upon 
phenomenologically motivated quark-antiquark potentials. There is empiri- 
cal evidence that this behavior persists for small loops, and even for the pla- 
quette, whose expectation value seems well described by perturbation theory 
times a factor of (1 — cra^), at least for /3 near 6.0 [Q. At (3 = 6.0, aa^ is 
about 0.04, making the nonperturbative part of the plaquette about 10% of 
the size of the 0{g'^) part. This is perhaps some indication of the relative 
sizes of perturbative and nonperturbative effects. 

To continue with our toy model, we assume that the mean-field parameter 
Mo inherits a nonperturbative contribution — from the plaquette: 

= „P-t - aaV4 . (101) 

This amounts to less than 2% of Uq^ aX (3 = 6.0, and less at higher /3's, so 
it has a negligible effect on renormalization constants like Zm and Z^. One 
place where one might worry about nonperturbative corrections is in the very 
divergent zero-point energy, 

Eq ^ (1-no) (l/a + S/Ma^) 

^ E^""'* + (Ta/4 + 3a/4M . (102) 

Here there is a contribution that does not vanish with a; factors of 1/a 
generated by perturbative power-law divergent loops cancel the factors of 
a in the nonperturbative term. However this contribution amounts to only 
about 30 MeV for b quarks, much less than the perturbative contribution of 
about 900 MeV at /5 = 6. 

Nonperturbative terms such as those in Eq. ( |102| ) for Eq mean that the 
perturbative relation, Eq. (pHD, between quark mass and meson mass cannot 
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be exact for any lattice spacing. But for (3 = 6.0, the uncertainty that resuhs 
is only a few percent of the b quark mass. Of course, ours is only a toy model 
for nonperturbative effects. We do not really know that the nonperturbative 
part of Mo is suppressed by rather than, say, or a^/^. But we do know 
that it is suppressed relative to the 0{g^) contributions. Thus we can use 
perturbative results to bound nonperturbative contributions; as long as per- 
turbative corrections are not individually large, nonperturbative corrections 
cannot be very important. In some cases, the use of the measured value of 
Mo, rather than the calculated one, may account for some of the nonpertur- 
bative physics. This is most likely the case, for example, when renormalizing 
cloverleaf E's and B's with the plaquette average. 

Finally it is worth noting that the nonperturbative contribution to the 
plaquette expectation value is almost comparable in magnitude to the O^g"^) 
perturbative contribution at /? = 6.0. Thus there seems little point in com- 
puting much beyond first or second order in g^. Precision beyond a few 
percent will probably require nonperturbative tuning of couplings. 

6 Conclusions 

Nonrelativistic QCD provides one of the most efficient frameworks for sim- 
ulating heavy quarks. While an approximation to QCD, it may be system- 
atically improved; this paper provides initial steps toward that goal. We de- 
veloped general power-counting rules which allowed us to assess the relative 
importance of various corrections, and which we can use to fine tune the the- 
ory for specific applications. We computed all of the leading order corrections 
required by relativity, both for spin-independent and spin-dependent inter- 
actions. The theory was then adapted to the lattice, including all leading- 
order corrections from finite lattice spacing. Finally, we presented a simple 
mean-field procedure that automatically incorporates the largest radiative 
corrections to the NRQCD coupling constants. Our mean-field analysis also 
allowed us to estimate the importance of nonperturbative contributions to 
the couplings. 

In NRQCD, heavy-quark propagators are determined by a simple evo- 
lution equation, avoiding the need for costly matrix inversion. Our fully 
corrected evolution equation was presented in Section ^73| . Using this equa- 
tion, together with the mean-field improvement procedure described in Sec- 
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tion 1^, should give results for ip and T mesons that are accurate to perhaps 
10% or better, depending upon the measurement. The largest remaining 
errors are probably due to uncalculated 0{as) corrections to the NRQCD 
couplings beyond the mean-field contributions, and to light-quark vacuum 
polarization. These errors can be removed in the near future. The coupling 
constant corrections are computed using ordinary one-loop weak-coupling 
perturbation theory; work has already begun on these. The dominant con- 
tribution from light-quark vacuum polarization is probably insensitive to 
light quark masses less than 100-200 MeV, and the extrapolation to realistic 
quark masses should be quite smooth. It is therefore feasible to include the 
light quarks with current lattice technology. Once these systematic effects 
have been removed, simulations should be accurate to a few percent, where 
nonperturbative contributions to the couplings may become important. 

Simulations using the techniques presented in this paper should, in the 
near future, produce accurate spectra, decay rates, wavefunctions, and other 
matrix elements for all of the lowest-lying mesons in the ip and T families. 
These techniques should also be useful in studies of D and B, as NRQCD 
propagators are more efficiently generated than relativistic propagators, and 
much less afflicted by noise than static propagators 0. With these methods 
we are entering a new era of high-precision lattice simulations of quantum 
chromodynamics. 
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Appendix 



This appendix discusses the construction in NRQCD of lattice operators 
which reproduce the corresponding continuum operators to any desired order 
in a, at tree level. We first describe a natural technique which works for 
all bilinears in fermion fields, but turns out to be inconvenient in actual 
simulations. Then we show how to circumvent this problem by starting with 
a convenient lattice operator and computing corrections to it. Finally, we 
discuss in some detail the corrections to the gluon action in the nonabelian 
case. 

Consider a general local, gauge covariant, continuum operator, bilinear 
in the fermion fields and with canonical dimension d = n + 3, 

Oix) = ^\x)Kni'E, B, D, A)V^(x) . (103) 

Since the field strength components are proportional to commutators of co- 
variant derivatives, such an operator is a linear combination of the tensors 

0,,,...,,M = ^Hx)D,, . . . D^Xx) . (104) 

It is straightforward to construct lattice expressions for these tensors, accu- 
rate to a specific order in a, provided one defines the link variables correctly 
in terms of the continuum gauge potentials, as in Eq. (|55|). The relationship 
between the continuum covariant derivatives and the corresponding finite 
differences on the lattice is then given by Eq. (0), and one can rewrite the 
tensors in Eq. ( p.04| ) using, for example, the symmetric expression 

^ 'ln(l + aA;+)) -ln(l-aA(-))l . (105) 



^ 2a 



In practice of course one always works with the first few terms in the power 
series expansion of Eq. (|105|). One can then use the fact that A|j+) = A|j~) + 



0{d) to construct an expression for that is maximally local on the lattice, 
as for example Eq. (^^ . To see what kind of lattice expressions are generated 
through this procedure, consider for example 

0[^y] {x) = - ig i)\x)F^^{x)i){x) . (106) 

Using Eq. (|57|), this becomes, up to corrections C(a^), 

0[h(x) = V^t(^)[(A(±)A(±) _^a(±)aWa(±)a(-) (107) 
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In terms of link variables the result is simpler and more transparent, and it 
is shown in Fig. 2. There one can see the main characteristic common to all 
lattice operators constructed this way: the quark field and its hermitian con- 
jugate are evaluated at different points on the lattice. In actual simulations, 
this turns out to be a serious disadvantage: it is much more economical to 
define both fermions at the same site. 

Fortunately, it is as easy to go from a lattice operator written in terms of 
link variables to the corresponding continuum expression. We can therefore 
start with any operator having the appropriate continuum limit, such as our 
cloverleaf F^l!{x), and then correct it to the desired order in a. Consider for 
example the first term in the cloverleaf definition of the field strength tensor, 
sandwiched between spinors: 

^Kx)U,,,U,^,f,,MU,,^,UlMx) . (108) 

It is easy to find an expression written in terms of lattice derivatives that 
will reproduce this term. One starts by replacing each matrix with a for- 
ward derivative A^j"^-*, and each matrix Uj^ with the corresponding backward 
derivative. This gives 

^t(x)AWAWA(-)A(-)V^(x) , (109) 

which contains Eq. ( |108| ), plus correction terms with at most three t/'s, that 
are easily calculated. The procedure is then iterated until all terms with a 
lower number of U's are canceled. Once the desired lattice operator is written 
in terms of lattice derivatives, the corresponding continuum operator can be 
calculated using, once again, Eq. (^), to the desired order in a. 

This procedure identifies the terms which need to be subtracted from the 
original definition. Once known in the continuum, it is easy to guess a lattice 
expression that will cancel them, yielding the corrected lattice operator. The 
result for the cloverleaf field strength is, as expected, given by Eq. ( |63D and 
Eq. (0). 

The correction to the gluon action in the nonabelian theory is somewhat 
more complicated to check directly. We start by noting that, to all orders in 
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a, the link variable can be written as 



x,ii 



exp 



CXJ 



.n=l 



(110) 



This is a consequence of the Baker-Campbell-Hausdorff theorem, that al- 
lows to rewrite the path-ordered exponential as exponential of a series of 
commutators. Our definition of the link variable, with the gauge potentials 
evaluated at mid-link, implies for example 



(111) 



u 



3,11 







< 


= = 




- ^^f,,,Al^,Ai 



where fbcd are the structure constants of SU{3). By the same token, we can 
write the plaquette matrix as 



Up(x,iiiy) = exp 



.n=l 



where, with the fields at the center, one easily finds that 



(112) 



(113) 



The series expansion in powers of a of the trace of Eq. (|112|) can be consid- 
erably simplified by using the trace identities 



Tr(Tj = 
Tr(T„T,Te: 



Tr(raTfc) = -5ah 

-;{dabc + ifabc) , 



(114) 



where dabc is the completely symmetric invariant tensor of SU{3). This yields 



P(x4ip) 



0,4 q5 
16 



(115) 



a 
24 



7\ 
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so that the only coefficients we need to calculate are Ca^^jy and C4,^jy. The 
calculation is rather lengthy, but can be performed using a symbolic manip- 
ulation program, such as "Mathematica" M. The result is 



[{d, + d,){A^ + A,),F^ 



+ -^[A^ + A^, [A^, + A^, F^^]] 



(116) 



Substituting Eq. ( |113| ) and Eq. ( |116| ) in Eq. ( |115| ) one can verify, as a nontriv- 
ial check, that all gauge non-invariant terms disappear, as do the coefficients 
of odd powers of a. The final result is, as expected 



Tr 



R (Up{x. 



^ TT{Fp - ^—Tt 



2 '"^ 24 
the nonabelian generalization of Eq. (I 



F,,iDl + Dl)F,, 



(117) 
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Operator 


Estimate 


Description 


V' 
X 
Dt 




quark (annihilation) field 
antiquark (creation) field 
gauge covariant time derivative 


D 


Mv 


gauge covariant spatial derivative 


g(t> 

gE 


Mv^ 
Mv^ 

M^v^ 


scalar potential (Coulomb gauge) 
vector potential (Coulomb gauge) 
chromoelectric field 


gB 




chromomagnetic field 



Table 1: The component fields and operators for the NRQCD action for 

heavy quarks. The estimated magnitude of each in a quarkonium state and 
in Coulomb gauge is given in terms of the quark mass M and typical velocity 

V. 
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aEo/g-": 



{Zm - I) 1 9 



2. 



Ma 


P.Th. 


M.F.Th. 


oo 


0.17 


0.08 


5 


0.21 


0.14 


2.5 


0.25 


0.18 


1.5 


0.30 


0.23 


oo 


0.07 


0.08 


5 


0.04 


0.03 


2.5 


0.07 


0.03 


1.5 


0.06 


0.0 


oo 


0.04 





5 


0.07 


0.05 


2.5 


0.10 


0.10 


1.5 


0.15 


0.17 



Table 2: The coefficient of in expansions of the NRQCD renormahzation 
parameters. Results coming from exact perturbation theory (P.Th.) are 
compared with estimates based on mean-field theory (M.F.Th.) for a variety 
of masses M. The parameters are appropriate to the evolution equation 
used by Davies and Thacker. For Zm, n = 1 for Ma — oo and 5; n = 2 for 
Ma — 2.5 and 1.5. Eq and are independent of n. An infrared divergence, 
which is shared with the continuum theory, has been omitted from Z^. 
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